%% calculate full lifetime earnings 

% as of child age 1 (t=1)
PDV_W_m=zeros(N,1);
PDV_W_f=zeros(N,1);

% as of child age 5 (t=5)
PDV_W_m2=zeros(N,1);
PDV_W_f2=zeros(N,1);

% as of child age 5 (t=5): only over ages 5-12
PDV_W_m3=zeros(N,1);
PDV_W_f3=zeros(N,1);

T_R=zeros(N,1);
T_D=zeros(N,1);

% 2002 wage relative to predicted wage
W_m_base=W_m./exp([expr_m expr_m.^2]*param_W_m);
W_f_base=W_f./exp([expr_f expr_f.^2]*param_W_f);

% 2002 age
age_base=age;

% experience at child birth
expr_m0=expr_m-age_base;
expr_f0=expr_f-age_base;

for i=1:N
    
    % average age of parents
    if (married(i)==0)
        age_p=age_m(i);
    else
        age_p=(age_m(i)+age_f(i))/2;
    end
    
    % year of retirement and death
    T_R(i)=J_R-(age_p-age(i));
    T_D(i)=J_D-(age_p-age(i));

    % loop over working years
    for t=1:T_R(i)

        x=expr_m0(i)+t;
        PDV_W_m(i)=PDV_W_m(i)+R^(-t+1)*W_m_base(i)*exp([x x^2]*param_W_m);
        
        if (t>=5)
            PDV_W_m2(i)=PDV_W_m2(i)+R^(-t+5)*W_m_base(i)*exp([x x^2]*param_W_m);
            if (t<=12)
                PDV_W_m3(i)=PDV_W_m3(i)+R^(-t+5)*W_m_base(i)*exp([x x^2]*param_W_m);
            end
        end
        
        if (married(i)==1)
            x=expr_f0(i)+t;
            PDV_W_f(i)=PDV_W_f(i)+R^(-t+1)*W_f_base(i)*exp([x x^2]*param_W_f);
            
            if (t>=5)
                PDV_W_f2(i)=PDV_W_f2(i)+R^(-t+5)*W_f_base(i)*exp([x x^2]*param_W_f);
                if (t<=12)
                    PDV_W_f3(i)=PDV_W_f3(i)+R^(-t+5)*W_f_base(i)*exp([x x^2]*param_W_f);
                end
            end
        end
    end
end

% update full income
W=(PDV_W_m+PDV_W_f)*5200;